Transcriptomics integrated with metabolomics reveals partial molecular mechanisms of nutritional risk and neurodevelopment in children with congenital heart disease

Purpose To explore molecular mechanisms affecting nutritional risk and neurodevelopment in children with congenital heart disease (CHD) by combining transcriptome and metabolome analysis. Methods A total of 26 blood and serum samples from 3 groups of children with CHD low nutritional risk combined with normal neurodevelopment (group A), low nutritional risk combined with neurodevelopmental disorders (group B) and high nutritional risk combined with normal neurodevelopment (group C) were analyzed by transcriptome and metabolomics to search for differentially expressed genes (DEGs) and metabolites (DEMs). Functional analysis was conducted for DEGs and DEMs. Further, the joint pathway analysis and correlation analysis of DEGs and DEMs were performed. Results A total of 362 and 1,351 DEGs were detected in group B and C compared to A, respectively. A total of 6 and 7 DEMs were detected in group B and C compared to A in positive mode, respectively. There were 39 and 31 DEMs in group B and C compared to A in negative mode. Transcriptomic analysis indicated that neurodevelopment may be regulated by some genes such as NSUN7, SLC6A8, CXCL1 and LCN8, nutritional risk may be regulated by SLC1A3 and LCN8. Metabolome analysis and joint pathway analysis showed that tryptophan metabolism, linoleic and metabolism and glycerophospholipid metabolism may be related to neurodevelopment, and glycerophospholipid metabolism pathway may be related to nutritional risk. Conclusion By integrating transcriptome and metabolome analyses, this study revealed key genes and metabolites associated with nutritional risk and neurodevelopment in children with CHD, as well as significantly altered pathways. It has important clinical translational significance.


Introduction
Congenital heart disease (CHD) is a congenital heart malformation caused by abnormal development of fetal heart and large blood vessels, and is the most common heart disease in children (1,2).CHD is the most common birth defect, and prevalence of CHD in live born infants is 1%-1.2%(1,3).Due to the variety, complexity and difficulty in treating CHD, its mortality is on the rise in China (4).CHD increases the risk of neurodevelopmental disorders, including cognitive, motor, social adjustment, and behavior disorders (5)(6)(7).Neurodevelopmental delays may reduce the lifespan of children with CHD and cause cognitive or intellectual impairment (8).This not only increase the risk of worse prognosis, but also has serious impact on the entire life, including academic performance, employment opportunities, psychological and overall quality of life (8,9).Hence, it is essential to explore the molecular mechanism of neurodevelopment in children with CHD and the interventions targeted at the early stage of the development to improve the quality of life.Growth failure is also a common problem in children with CHD (10,11).The etiology and pathological mechanism of slow physical and mental development in infants with CHD remain unclear, but there are many risk factors associated with it, including severity of cyanosis, chronic hypoxia, hemodynamic changes, repeated infections and heart failure, repeated hospitalization, and feeding disorders (11,12).Nutrition, parents' education level, living environment, social and family environment, living style and economic level will also affect the physical and intellectual development of children with CHD (11,13).Perioperative young children required more caloric and nutrient intake to promote adequate growth and psychomotor development, and energy and protein deficiencies can increase infections and inflammatory responses, impair wound healing, prolong hospital stays, and may increase the incidence of postoperative complications (14,15).Therefore, optimal nutrition is crucial to improve the short-and long-term prognosis of children with CHD.
Omics analysis techniques are increasingly used to identify potential biomarkers and elucidate causes and mechanisms associated with disease (16,17).Metabolomics can analysis small molecule metabolites to reflect the biological metabolic characteristics of disease states and is used extensively in biomarker discovery (18).Dong S et al. analyzed the heart's metabolic remodeling to hypoxia using metabolomics and found protein synthesis and aerobic energy production were reduced in patients with cyanotic CHD and NAD may play an important role in response to hypoxia (19).Jin N et al. found S-Adenosyl methionine, guanine and N-terminal pro-brain natriuretic peptide could be used as biomarkers for pulmonary arterial hypertension associated with CHD from CHD (20).Transcriptome research can study gene function and gene structure at the whole level and reveal the molecular mechanism of specific biological processes and disease occurrence (21).Liu G et al. identified genes and enriched pathways associated with radiotherapy in nasopharyngeal carcinoma, revealing the molecular mechanism of radiotherapy resistance, which is helpful for future research on radiotherapy resistance function (22).The combined transcriptomics and metabolomic analysis of high-risk neuroblastoma and low and intermediate-risk neuroblastoma identified 4 aberrant pathways and developed a risk classification diagnostic model, providing insights for high-risk neuroblastoma early diagnosis (23).Therefore, transcriptomics and metabolomics play an important role in diagnosis of diseases, and their molecular mechanisms can be deeply explored.
In this study, transcriptomic and metabolomic analysis of blood and serum samples from 26 children with CHD was performed to search for key DEGs and DEMs, and explore molecular mechanisms affecting nutritional risk and neurodevelopment in children with CHD.
2 Methods and materials

Sample collection
A total of 26 blood and serum samples were collected from hospitalized children with CHD in the Dalian Children's Hospital.Blood samples (2 ml/sample) were used for transcriptomic analysis and serum samples (1 ml/sample) were used for metabolomics analysis.Main inclusion criteria: (1) confirmed diagnosis of CHD; (2) aged 0-24 months; (3) written informed consent was obtained from child's guardians.Main exclusion criteria: (1) children with nutritional impairments due to major non-cardiac diseases; (2) multiple pregnancy; (3) chromosomal abnormality; (4) structural brain malformation; (5) placental dysfunction; (6) intrauterine growth retardation of the fetus; (7) received the prescribed nutrition supply and blood transfusion half a month before hospitalization.Blood and serum samples were collected 24 h before discharge from the hospital (or after surgery).The children with CHD were divided into three groups according to the nutritional risk and neurodevelopmental status: low nutritional risk combined with normal neurodevelopment (n = 14, group A), low nutritional risk combined with neurodevelopmental disorders (n = 5, group B), high nutritional risk combined with normal neurodevelopment (n = 7, group C).The screen tool for risk on nutritional status and growth (STRONGkids) was used for nutritional risk screening and nutritional assessment of hospitalized children over 1 month.The Revised Gesell Developmental Schedules were used to assess the neurodevelopmental status of CHD children aged 6 months to 24 months.This study was approved by the Ethics Committee of Dalian Children's Hospital (19015).

Transcriptomics detection through RNA-Sequencing (RNA-Seq) analysis
Total RNA was extracted from the samples with TRIzol.Sequencing By Synthesis (SBS) technology was used to sequence the cDNA library using Illumina's high-throughput sequencing platform.The raw data obtained by sequencing were converted into FASTQ format sequence data by base calling.The raw data was analyzed using a bioinformatic workflow, which includes quality control with MultiQC, adapter and quality trimming with Cutadapt, transcript quantification with featureCounts (GRCh38.primary_assembly.genome.fa).Differential expression analysis was performed using the DESeq2 in R package (4.0.5) (p value < 0.05 and |log 2 foldchange| > 1).The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) function enrichment analysis of differentially expressed genes (DEGs) were performed using the David database (https://david.ncifcrf.gov/tools.jsp).

Non-targeted metabolomics analysis via liquid chromatograph with tandem mass spectrometer (Lc-Ms/Ms)
The 300 μl acetonitrile was added to 100 μl serum at 4°C.After 20 min of centrifugation at 12,000 r/min (4°C), 100 μl supernatant was carefully extracted.Take 50 μl as quality control (QC) samples.All of the extracts were analyzed using LC-MS/MS technology with a high performance liquid chromatography (HPLC) and a high resolution mass spectrometer (HRMS).The chromatographic separation was performed on a Waters HSS T3 column.Metabolomics data were collected as follows: a column temperature of 45°C and the flow rate of 0.3 ml/min.mobile phase: A = 0.1% formic acid water, B = 0.1% formic acid acetonitrile.Mass spectrometry (MS) negative and positive mode conditions: sheath gas flow rates 30 arb, atomization pressure 50 psi.capillary voltage: positive mode 5,500 V, negative mode −4,500 V. Based on the MS detection, the original files were imported into Progenesis QI (Waters) software for data pre-processing and identification, and then the data quality control analysis was carried out to ensure the accuracy and reliability of the data.Multivariate statistical analysis was performed on the data to reveal the differences of metabolites of different components, the biological significance of metabolites was discovered through functional analysis.Orthogonal partial least squares discrimination analysis (OPLS-DA) was used to filter signals that are not relevant to classification.Potential metabolites were analyzed according to predicted value (VIP) and significance of variables in Student's t-test.VIP > 1 and p < 0.05 were considered as statistically significant.Differentially expressed metabolites (DEMs) were selected using KEGG (https://www.genome.jp/kegg/)and human metabolome database (HMDB) (https://hmdb.ca/), and further annotated in the KEGG com-pound database.The annotated metabolic pathways were classified according to the KEGG pathway database (https://www.kegg.jp/kegg/pathway.html).

Joint analysis of the metabolomics and transcriptomics
Combined pathway analysis of DEMs and DEGs were performed using MetaboAnalyst 6.0 (https://www.metaboanalyst.ca/).The cor. test function in R software was used to calculate Pearson correlation coefficient to construct the correlation analysis between the above metabolites and genes.

Transcriptome sequencing analysis and identification of DEGs
To explore the correlation between nutritional risk and neurodevelopment in children with CHD, we performed transcriptomic analysis of the three groups.Principal component analysis (PCA) was used to test the reliability of the experiment and the rationality of sample selection.Pearson correlation coefficient, R 2 > 0.8 was considered as good repeatability sample.As shown in Supplementary Figure S1, samples within the group have good repeatability, and samples between the groups have large differences.Sample N14 was significantly different from other samples, so it was excluded in the subsequent analysis.
A total of 362 DEGs were screened in group B compared to A, of which 145 up-regulated genes and 217 down-regulated genes were detected (Figures 1A,B).The top 20 up-regulated and down-regulated genes in group B compared to A are shown in Table 1.A total of 1,351 DEGs were screened in group C compared to A, of which 1,079 up-regulated genes and 272 down-regulated genes were detected (Figures 1C,D).The top 20 up-regulated and down-regulated genes in group C compared to A are shown in Table 2.

GO analysis and KEGG pathway analysis
GO and KEGG enrichment analysis of the obtained DEGs was performed to analyze the key functions and metabolic pathways (Figure 2; Supplementary Table S1).The top three significantly enriched terms were innate immune response, response to virus and defense response to virus in biological processes in group B compared with group A. The top three significantly enriched terms were extracellular region, extracellular space and plasma membrane in cellular component in group B compared with group A. The top three significantly enriched terms were laminin binding, oligoadenylate synthetase activity and CXCR chemokine receptor binding in molecular function in group B compared with group A (Figure 2A).The top three significantly enriched terms were innate immune response, defense response to virus and inflammatory response in biological processes in group C compared with group A. The top three significantly enriched terms were plasma membrane, membrane and tertiary granule membrane in group C compared with group A. The top three significantly enriched terms were protein binding, NAD + nucleosidase activity and NAD + nucleotidase, cyclic ADP-ribose generating in molecular function in group C compared with group A (Figure 2C).KEGG analysis found that the top three altered pathways were influenza A, viral protein interaction with cytokine and cytokine receptor and cytokine-cytokine receptor interaction in group B compared with group A (Figure 2B), and the top three altered pathways were NOD-like receptor signaling pathway, cytokine-cytokine receptor interaction and osteoclast differentiation in group C compared with group A (Figure 2D).In summary, compared with group A, the key functions of both groups B and C are significantly enriched innate immune response and defense response to virus and plasma membrane, and the significant enrichment pathway is cytokine-cytokine receptor interaction.Metabolomics analysis and identification of DEMs.
Correlation analysis of the QC samples showed that R 2 values were greater than 0.8 (Supplementary Figures 2A,B), indicating that the entire analysis process was stable and reproducible.OPLS-DA was performed to screen the reliable metabolites that lead to the classification difference.The results showed that the metabolites of the two groups with low nutritional risk combined with neurodevelopmental disorders and high nutritional risk combined with normal neurodevelopment were significantly different from those of the healthy control group in both positive and negative modes (Supplementary Figures 2C-F).
Metabolites with VIP > 1 and p < 0.05 in QC samples were selected as DEMs.In positive mode, 6 DEMs were identified in group B compared to A, all of which were up-regulated metabolites (Figure 3A; Table 3).A total of 39 metabolites changed significantly in negative mode in group B compared to A, including 8 up-regulated and 31 down-regulated metabolites (Figure 3B; Table 4).Similarly, compared with A, we found a total of 7 DEMs in positive mode in group C, 4 of which were up-regulated and 3 down-regulated (Figure 3C; Table 3).In   3D; Table 4).

KEGG pathway analysis of DEMs
The DEMs were subjected to KEGG-enrichment analysis, and enriched pathways with significant differences were identified.The enriched metabolic pathways were Parkinson's disease, linoleic and metabolism and glycerophospholipid metabolism in group B compared to A (Figure 4A).The enriched metabolic pathway was glycerophospholipid metabolism in group C compared to A (Figure 4B).Compared with group A, the common enrichment metabolic pathway of group B and C was glycerophospholipid metabolism.Correlation analyses between transcriptomics and metabolomics.
To explore the factors affecting nutritional risk and neurodevelopment in children with CHD by linking important metabolites and genes through shared metabolic pathways, joint pathway analysis between DEGs and DEMs were performed.The pathways are shown in Figures 5A,B, which includes tryptophan metabolism and mucin type O-glycan biosynthesis in group B compared to A. Nine significantly altered pathways were revealed in group C compared to A, include glycerophospholipid metabolism, neomycin, kanamycin and gentamicin biosynthesis, nitrogen metabolism, glycerolipid metabolism, starch and sucrose metabolism, arginine biosynthesis, mucin type O-glycan biosynthesis, glycosaminoglycan biosynthesis-heparan sulfate/ heparin and galactose metabolism (Figures 5C,D).Correlation analysis utilized Pearson calculation to show the correlation of the DEMs and DEGs. Figure 6 showed strong correlations among transcripts and metabolites.

Discussion
In this study, we divided children with CHD into three groups based on nutritional risk and neurodevelopment, and integrated transcriptomics and metabolomics data to explore molecular mechanisms affecting nutritional risk and neurodevelopment in children with CHD.We found that NSUN7, SLC6A8, CXCL1 and LCN8 were significantly different in groups B compared to A, and SLC1A3 and LCN8 were significantly different in groups C compared to A. NSUN7, a m 5 C RNA methyltransferase gene, is differentially expressed in Alzheimer's disease patients (24).Solute carrier (SLC) are a group of membrane transport proteins that facilitate the transport of various substrates across biofilms (25).SLC6A8 could modulate human creatine levels and suppress colon cancer progression (26).Previous studies have shown that SLC6A8 can be up-regulated by p65/NF-κB transcription and mediate intracellular creatine accumulation in hypoxia (27).Creatine deficiency may be manifested as intellectual and behavioral disorders (28), so we suspect that SLC6A8 is related to neurodevelopment in children with CHD.SLC1A3 is an aspartate/glutamate transporter that maintains electron transport chain and tricarboxylic acid cycle activity (29).Studies have also shown that SLC1A3 may affect glucose metabolism by activating the PI3 K/AKT signaling pathway to provide energy for cell proliferation (30), which also seems to confirm its correlation with nutritional risk.Besides, SLC1A3 was also related to the development of stress and depression (31).CXCL1 was expressed in nervous systems and involved in the development of inflammation and pain (32).LCN8, a member of the liposome protein family, was found to be associated with the severity of post-traumatic stress disorder symptoms by genomewide DNA methylation analysis (33).Therefore, we deduced that NSUN7, CXCL1, SLC6A8 and LCN8 may be involved in the neurodevelopment in children with CHD, SLC1A3 and LCN8 may be involved in the nutritional risk.We also found that cytokine receptor interaction was found to be the most altered pathway in both two groups compared to the control group.
Cytokines play an important regulatory role in various processes, including immune function, inflammation, hematopoiesis, cell growth and differentiation (34).Cytokine receptor interaction provides a new perspective to explore the molecular mechanism affecting nutritional risk and neurodevelopment in children with CHD.Metabolomics analysis found the levels of Albanol B, dihomolinoleic acid, lysophosphatidylethanolamine (LysoPE) and lysophosphatidylcholine (LysoPC) were significantly altered in group B compared to A, and phospholipids and amino acids were significantly altered in group C compared to A. Albanol B, a compound isolated in root bark, could be involved in Anti-Alzheimer's disease (35).Dihomo-linoleic acid could drive ferroptosis-mediated neurodegeneration (36).LysoPE and LysoPC are lysophospholipids, a type of phospholipid.Lysophospholipids are involved in many physiological processes, including regulating neuropathic pain, disorders of neuroectodermal development, and disorders of neuroectodermal development (37,38).Lipids are closely related to Alzheimer's disease (39,40).Several key pathways were identified in group B compared to A (neurodevelopment disorder enriched comparison), including Parkinson's disease, linoleic acid metabolism and glycerophospholipid metabolism.Glycerophospholipid metabolism is also a critical pathway in group C compared to A. Linoleic acid is essential for optimal growth and brain development in infants and can regulate key neurodevelopmental processes (41).Studies have shown that cancer patients will suffer from malnutrition, and the increase of nutritional risk will increase the content of linoleic acid in fecal supernatant, indicating that linoleic acid is a biomarker of high nutritional risk (42).Linoleic acid can synthesize polyunsaturated fatty acids required by human brain cells, which is conducive to promoting brain development and plays a certain role in promoting the development of children's intelligence (41).This is   consistent with our study that the linoleic acid metabolic pathway is associated with neurodevelopment in children with CHD.A basic study has demonstrated that dysregulation of glycerophospholipid and sphingolipid metabolism may have a negative impact on neurodevelopment in offspring, while lipid changes may disturb phospholipid homeostasis, altering membrane integrity, orientation, permeability, and function, leading to neurological dysfunction and degeneration (43).Choline is the main component of synthetic phospholipids, phosphatidylcholine and sphingomyelin, and plays an important role in neurogenesis and neural migration during fetal development (44).Gut microbiota plays an important role in the onset of depression, and metabolomics analysis showed that significant differences in glycerophospholipids and fatty acids metabolism between depressed mice and normal mice (45).Phospholipid metabolism plays an important role in energy metabolism throughout the body, and phospholipid synthesis is essential for normal development and health (46).It can be concluded that linoleic acid metabolism and glycerophospholipid metabolism may be involved in the nutritional risk and neurodevelopment of children with CHD.
In joint analysis, glycerophospholipid metabolism was found to be a key pathway in group C compared to A (nutritional risk enriched comparison), which was consistent with metabolomic results.In the network diagram of glycerophospholipid  Frontiers in Cardiovascular Medicine metabolism, the expression of GPD2 and AChE were up-regulation in group C compared to A, which was consistent with transcriptomic results.GPD2, a glycerol 3-phosphate dehydrogenase, was a component of the glycerol phosphate shuttle process that promotes glucose oxidation, the production of acetyl coenzyme A, the acetylation of histones, and the induction of genes encoding inflammatory mediators (47).AChE is a serine protease that hydrolyzes the neurotransmitter acetylcholine into acetate and choline, thus terminating neurotransmission (48).Tryptophan metabolism was found to be a key pathway in group B compared to A (neurodevelopment disorder enriched comparison).Tryptophan can be involved in the pathophysiology of different neuropsychiatric diseases through the serotonin pathway to produce serotonin as a neurotransmitter and melatonin as a neuromodulator (49).
Tryptophan metabolism is directly or indirectly regulated by intestinal microorganisms, and its metabolites have immune, metabolic and neuroregulatory functions, which has become a therapeutic target for various diseases (50).Tryptophan metabolism pathways are also associated with symptoms and neurodevelopment in children with autism spectrum disorders (51), which was consistent with our results.In summary, the content of metabolites and gene expression in blood and serum samples were analyzed to explore molecular mechanisms affecting nutritional risk and neurodevelopment in children with CHD.We further identified NSUN7, SLC6A8, CXCL1 and LCN8 as possible molecular markers associated with neurodevelopment and SLC1A3 and LCN8 were associated with nutritional risk.Linoleic acid metabolism, glycerophospholipid metabolism and tryptophan metabolism pathways were associated with neurodevelopment and glycerophospholipid metabolism pathway was associated with nutritional risk in children with CHD.By integration of metabolomics and transcriptomics data, our study provided the molecular basis for nutritional risk and neurodevelopment in children with CHD that could provide reference for clinical decision-making for individualized nutrition supply and demand in Chinese patients with CHD.This study also lays the groundwork for validation studies in longitudinal assessments of neurodevelopmental outcomes and preclinical investigations of therapeutic interventions targeting the identified molecular pathways.These findings need to be validated in large numbers of patients in the future and additional biological validation of the particular molecular processes in in vitro or in vivo investigations is necessary.

FIGURE 1
FIGURE 1 Identification of DEGs.(A) The volcano plot shows DEGs between group B and A. (B) The heatmap shows DEGs between group B and A. (C) The volcano plot shows DEGs between group C and A. (D) The heatmap shows DEGs between group C and A. DEGs, differentially expressed genes.

FIGURE 2
FIGURE 2 Go and KEGG analysis of DEGs.(A) GO analysis of biological processes, cellular components and molecular functions of DEGs between group B and A. (B) Pathway enrichment of the DEGs between group B and A. (C) GO analysis between group C and A. (D) Pathway enrichment of the DEGs between group C and A. The dot color represents the p-value and the dot size represents the number of DEGs.DEGs, differentially expressed genes.

FIGURE 3
FIGURE 3 Identification of DEMs.(A) The heatmap shows DEMs in positive mode between group B and A. (B) The heatmap shows DEMs in negative mode between group B and A. (C) The heatmap shows DEMs in positive mode between group C and A. (D) The heatmap shows DEMs in negative mode between group C and A. DEMs, differentially expressed metabolites.

FIGURE 4 KEGG
FIGURE 4 KEGG analysis in metabolomics.(A) Pathway enrichment of the DEMs between group B and A. (B) Pathway enrichment of the DEMs between group C and A. The dot color represents the p-value and the dot size represents the number of DEMs.DEMs, differentially expressed metabolites.

FIGURE 5 Joint
FIGURE 5 Joint-pathway analysis of DEGs and DEMs.(A) Joint-pathway analysis of DEGs and DEMs between group B and A. (B) Tryptophan metabolism pathway.(C) Joint-pathway analysis of DEGs and DEMs between group C and A. (D) Glycerophospholipid metabolism pathway.DEGs, differentially expressed genes; DEMs, differentially expressed metabolites.

FIGURE 6
FIGURE 6 Correlation analysis of DEGs and DEMs between group B and A (A), group C and A (B) the quadrangle indicates metabolites, and the nodes indicate genes.

TABLE 1
The top 20 up-regulated and down-regulated genes in group B compared to A.

TABLE 3
Differentially expressed metabolites in B vs. A and C vs. A in positive mode.

TABLE 4
Differentially expressed metabolites in B vs. A and C vs. A in negative mode.